{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "SgQuIgMpD084",
      "metadata": {
        "id": "SgQuIgMpD084"
      },
      "source": [
        "# Cell type identification in scRNA-seq data using ACTINN and Deepchem"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "Q8XrmxgyLkQu",
      "metadata": {
        "id": "Q8XrmxgyLkQu"
      },
      "source": [
        "## Why is RNA important?\n",
        "\n",
        "To understand what scRNA-seq(Single-cell RNA Sequencing) data represents, let's first explore the role of RNA (Ribonucleic Acid) inside cells. All cells in our body, regardless of whether they make up tissues in the heart, brain, or skin, contain the same genetic code (DNA). Yet, despite having identical DNA, cells look different, perform distinct functions, and form specialized tissues and organs.\n",
        "\n",
        "How is this possible?\n",
        "\n",
        "The key lies in differential [gene expression](https://en.wikipedia.org/wiki/Gene_expression), the process by which genetic instructions are used to synthesize gene products like proteins. DNA acts as the blueprint, but it is RNA, specifically messenger RNA (mRNA), that translates these instructions into action. Different types of cells transcribe distinct subsets of genes into mRNAs, leading to the synthesis of unique sets of proteins. These distinct protein profiles allow cells to specialize and perform various functions within the body [1,2].\n",
        "\n",
        "Messenger RNA (mRNA) acts as a critical intermediary, carrying genetic instructions from DNA to ribosomes,the cellular machinery responsible for protein synthesis.\n",
        "\n",
        "Proteins like enzymes drive metabolic pathways, while structural proteins such as actin and myosin give muscles their ability to contract and move. Cytoskeletal proteins help maintain cell shape and integrity.\n",
        "\n",
        "Therefore, the unique set of mRNA molecules present in a cell at any given moment(its **transcriptome**) reveals which genes are \"turned on.\" This transcriptome dictates which proteins the cell is making, which in turn defines the cell's type and function. Measuring this transcriptome is the fundamental goal of single-cell RNA sequencing (scRNA-seq).\n",
        "\n",
        "A high-level overview of single-cell RNA sequencing (scRNA-seq) is illustrated in the image below.\n",
        "\n",
        "<img src=\"https://learn.gencore.bio.nyu.edu/wp-content/uploads/2018/01/scRNA-overview.jpg\" width=800 height=600>\n",
        "\n",
        "For a more in-depth understanding of the various techniques used in scRNA-seq and the challenges associated with the field, refer to the comprehensive review by Hwang et al. (2018) published in Experimental & Molecular Medicine.[5]\n",
        "\n",
        "## Understanding RNA-seq Platforms\n",
        "RNA sequencing is a method used to measure which genes are active in a cell or tissue. Depending on the question you're trying to answer, you might use bulk RNA-seq, which measures average gene activity across many cells, or single-cell RNA-seq, which measures gene expression in individual cells. A RNA-seq platform is the technology or system used to perform RNA sequencing, the process of reading and measuring the RNA molecules present in a biological sample.\n",
        "\n",
        "\n",
        "Bulk RNA-seq Platforms\n",
        "Bulk RNA-seq gives you a snapshot of average gene expression from a population of cells. It’s commonly used to compare conditions like healthy vs. diseased tissue or treated vs. untreated samples.\n",
        "\n",
        "- Illumina\n",
        "\n",
        "  - Illumina is the most widely used platform for bulk RNA-seq.\n",
        "  - It generates short-read sequences (usually 50–150 base pairs) with very high accuracy.\n",
        "  - It uses a method called sequencing by synthesis, where fluorescent signals are recorded as DNA strands are built one base at a time.\n",
        "  - This platform is ideal for large-scale experiments, and it's supported by many bioinformatics tools.\n",
        "\n",
        "\n",
        "Single-cell RNA-seq Platforms\n",
        "Single-cell RNA-seq lets researchers study gene expression in individual cells, which is essential for identifying rare cell types, understanding cellular diversity, and mapping developmental or disease processes.\n",
        "\n",
        "- 10x Genomics Chromium\n",
        "\n",
        "  - 10x Genomics is the leading commercial platform for single-cell RNA-seq.\n",
        "  - It uses droplet microfluidics, where each cell is captured in a tiny droplet with a barcode that labels its RNA.\n",
        "  - This allows the analysis of thousands to tens of thousands of cells at once.\n",
        "  - It’s ideal for studying complex tissues like tumors, brains, or blood, where many different cell types are present.\n",
        "\n",
        "- SMART-Seq\n",
        "\n",
        "  - SMART-Seq is a more sensitive method that captures full-length transcripts from single cells.\n",
        "  - Instead of droplets, it usually uses well plates to handle cells.\n",
        "  - It has lower throughput compared to 10x (hundreds vs. thousands of cells) but gives richer information per cell.\n",
        "  - This is useful when you're studying small, precise populations, such as stem cells or neurons, and need full transcript coverage.[3]\n",
        "\n",
        "## Libraries commonly used for scRNA-seq analysis\n",
        "\n",
        "[ScanPy](https://scanpy.readthedocs.io/en/stable/) is a scalable toolkit for analyzing single-cell gene expression data. It includes methods for preprocessing, visualization, clustering, pseudotime and trajectory inference, differential expression testing, and simulation of gene regulatory networks.\n",
        "\n",
        "[ANNDATA](https://github.com/theislab/anndata) was presented alongside ScanPy as a generic class for handling annotated data matrices that can deal with the sparsity inherent in gene expression data.\n",
        "\n",
        "This is also a good point to explore DeepChem's [tutorial](https://deepchem.io/tutorials/scanpy/) on building an scRNA-seq analysis pipeline using ScanPy.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8e4b19fb",
      "metadata": {
        "id": "8e4b19fb"
      },
      "source": [
        "## ACTINN (automated identification of cell types in single cell RNA sequencing )\n",
        "\n",
        "This tutorial guides you through working with single-cell RNA sequencing (scRNA-seq) data and demonstrates how to train and evaluate a model using Deepchem. The goal is to replicate the cell type identification experiments presented in the ACTINN.[4]\n",
        "\n",
        "**Cell Type Identification**\n",
        "\n",
        "While scRNA-seq reveals which genes are active in each individual cell, it doesn't tell us what type of cell it is. Identifying cell types based on gene expression profiles unlocks a wide range of biological insights and enables powerful downstream analyses. Below are some key applications:\n",
        "\n",
        "- **Understanding tissue composition:** Tissues and organs are made up of multiple cell types working together.scRNA-seq allows researchers to \"deconstruct\" a tissue into its cellular building blocks and determine the proportions of each cell type.\n",
        "- **Comparing Conditions (e.g., Healthy vs. Diseased):** Understanding differences between conditions often requires comparing the same cell types across them. Cell type identification lets researchers perform differential expression or pathway analysis within each cell type.Example: Comparing gene expression in immune cells from healthy vs. cancerous tissue to find therapeutic targets.\n",
        "- **Characterizing the Tumor Microenvironment (TME)**\n",
        "Tumors are not just made of cancer cells, they’re surrounded by immune cells, blood vessels, and connective tissue. Cell type identification helps map out this tumor microenvironment, showing how each cell type contributes to tumor growth or suppression.\n",
        "- **Predicting and Improving Immunotherapy Response**\n",
        "Immunotherapies (like checkpoint inhibitors) work by reactivating the immune system to fight cancer. Cell type identification helps predict which patients are more likely to respond to these treatments, based on the immune cell composition of their tumors. Example: Tumors with a higher abundance of active CD8+ T cells often respond better to certain immunotherapies.[9]\n",
        "\n",
        "\n",
        "**Traditional approach:**\n",
        "\n",
        "To identify cell types from scRNA-seq data, typically cells are first grouped into different clusters in an unsupervised way, and the number of clusters allows us to approximately determine how many distinct cell types are present in the sample. Each cluster should contain cells with similar expression profiles, and so the aggregated profile of a cluster increases the signal to noise of the expression estimates. To attempt to interpret the identity of each cluster, marker genes are found as those that are uniquely highly expressed in a cluster, compared to all the other clusters. These canonical markers are then used to assign the cell types for the clusters, by cross referencing the markers with lists of previously characterized cell type specific markers. While this process is able to identify cell types, there are some limitations: \n",
        "\n",
        "**Limitations:**\n",
        "1.  Since the clustering method is unsupervised, all sources of variation influence the formation clusters, including effects that are not directly related to cell types such as differential expression induced by cell cycles. \n",
        "2. It is often difficult to find an optimal match between the marker genes associated with each cluster and the canonical markers( well established markers used to identify and classify different cell types) for specific cell types. Moreover, depending on the clustering parameters used, one cluster might contain multiple cell types, or one cell type could be split into multiple clusters. \n",
        "3. Using canonical markers to assign cell types requires background knowledge of cell type specific markers, and sometimes these are not well characterized or difficult to find in the literature. Moreover, some canonical markers may be expressed by more than one cell type, and some cell types may have no known markers. \n",
        "4. The same types of cells processed by two distinct scRNA-seq techniques tend to cluster separately due to technical batch effects, which complicates cell type identification in composite datasets. \n",
        "5. Cell subtypes are often very similar to each other, which limits efforts to separate them accurately into different clusters.\n",
        "\n",
        "To overcome many of the limitations of traditional cell type identification approaches described above, **ACTINN** uses a neural network with three hidden layers, trained it on scRNA-seq datasets with predefined cell types, and predicted cell types in other datasets based on the trained parameters.The input layer has a number of nodes equal to the number of genes in the training set. The 3 hidden layers have 100, 50 and 25 nodes, respectively. The output layer has a number of nodes equal to the number of cell types in the training set.[4]\n",
        "\n",
        "**Machine Learning based approach:**\n",
        "- Machine learning is a way for computers to learn patterns from data without being explicitly programmed. In the case of scRNA-seq, ML models can learn how gene expression patterns correspond to specific cell types, just like how facial recognition software learns to identify people from photos.\n",
        "- Instead of relying on a few known marker genes, ML models learn from thousands of genes at once, detecting subtle combinations that define each cell type. Once trained, the model can label thousands or millions of cells in seconds, making it ideal for large-scale experiments.\n",
        "\n",
        "\n",
        "## Colab\n",
        "\n",
        "This tutorial and the rest in this sequence can be done in Google colab. If you'd like to open this notebook in colab, you can use the following link.\n",
        "\n",
        "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/deepchem/deepchem/blob/master/examples/tutorials/Cell_type_identification_using_scRNAseq_data.ipynb)\n",
        "\n",
        "## Setup\n",
        "\n",
        "To run DeepChem within Colab, you'll need to run the following installation commands. You can of course run this tutorial locally if you prefer. In that case, don't run these cells since they will download and install DeepChem in your local machine again."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "587b693f",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 766
        },
        "id": "587b693f",
        "outputId": "ee7a6e3e-e830-43eb-ab20-ffd7d3dd9fb0"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "No normalization for SPS. Feature removed!\n",
            "No normalization for AvgIpc. Feature removed!\n",
            "No normalization for NumAmideBonds. Feature removed!\n",
            "No normalization for NumAtomStereoCenters. Feature removed!\n",
            "No normalization for NumBridgeheadAtoms. Feature removed!\n",
            "No normalization for NumHeterocycles. Feature removed!\n",
            "No normalization for NumSpiroAtoms. Feature removed!\n",
            "No normalization for NumUnspecifiedAtomStereoCenters. Feature removed!\n",
            "No normalization for Phi. Feature removed!\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /opt/anaconda3/envs/dc_new/lib/python3.10/site-packages/tensorflow/python/util/deprecation.py:588: calling function (from tensorflow.python.eager.polymorphic_function.polymorphic_function) with experimental_relax_shapes is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "experimental_relax_shapes is deprecated, use reduce_retracing instead\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "Skipped loading modules with pytorch-geometric dependency, missing a dependency. No module named 'torch_geometric'\n",
            "Skipped loading modules with pytorch-geometric dependency, missing a dependency. cannot import name 'DMPNN' from 'deepchem.models.torch_models' (/Users/harin/Desktop/deepchem/deepchem/deepchem/models/torch_models/__init__.py)\n",
            "Skipped loading modules with pytorch-lightning dependency, missing a dependency. No module named 'lightning'\n",
            "Skipped loading some Jax models, missing a dependency. No module named 'jax'\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "'2.8.1.dev'"
            ]
          },
          "execution_count": 1,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "!pip install deepchem\n",
        "import deepchem as dc\n",
        "dc.__version__"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fd13efd5",
      "metadata": {
        "id": "fd13efd5"
      },
      "source": [
        "## Data Loading"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ec53f3e8",
      "metadata": {
        "id": "ec53f3e8"
      },
      "source": [
        "A subset of the 'tma_ss2_cleaned' dataset from this source[8], which was originally used in ACTINN, is utilized here.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ed9a6785",
      "metadata": {},
      "source": [
        "**Data Splitting**\n",
        "\n",
        "There are two approaches to creating the train and test sets:\n",
        "1. Splitting a single dataset into training and testing portions.\n",
        "2. Using datasets from different sources as the train and test sets."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "eccc26d8",
      "metadata": {
        "id": "eccc26d8"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "current_dir = os.path.dirname(os.path.realpath('__file__'))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "1539bd90",
      "metadata": {
        "id": "1539bd90"
      },
      "outputs": [],
      "source": [
        "file_names = ['train_sample.h5', 'train_label_sample.txt', 'test_sample.h5', 'test_label_sample.txt']\n",
        "base_url = 'https://raw.githubusercontent.com/deepchem/deepchem/master/examples/tutorials/assets/scRNAseq'\n",
        "\n",
        "for i in range(len(file_names)):\n",
        "  dc.utils.download_url(\n",
        "      f\"{base_url}/{file_names[i]}\",\n",
        "      current_dir,\n",
        "      f\"{file_names[i]}\"\n",
        "  )"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "25531170",
      "metadata": {
        "id": "25531170"
      },
      "source": [
        "### 1. Same source"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "797fca67",
      "metadata": {
        "id": "797fca67"
      },
      "outputs": [],
      "source": [
        "labels = pd.read_csv(os.path.join(current_dir,'train_label_sample.txt'), header=None, sep='\\t') \n",
        "dataset = pd.read_hdf(os.path.join(current_dir,'train_sample.h5'))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "GyJFjnTixWfW",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 443
        },
        "id": "GyJFjnTixWfW",
        "outputId": "5c45f9ec-8f5f-4a70-c5d8-ee2a59053456"
      },
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>tma_facs_14301</th>\n",
              "      <th>tma_facs_4756</th>\n",
              "      <th>tma_facs_14382</th>\n",
              "      <th>tma_facs_423</th>\n",
              "      <th>tma_facs_15612</th>\n",
              "      <th>tma_facs_4729</th>\n",
              "      <th>tma_facs_8174</th>\n",
              "      <th>tma_facs_9855</th>\n",
              "      <th>tma_facs_9055</th>\n",
              "      <th>tma_facs_8668</th>\n",
              "      <th>...</th>\n",
              "      <th>tma_facs_5182</th>\n",
              "      <th>tma_facs_19106</th>\n",
              "      <th>tma_facs_11386</th>\n",
              "      <th>tma_facs_14006</th>\n",
              "      <th>tma_facs_7477</th>\n",
              "      <th>tma_facs_15352</th>\n",
              "      <th>tma_facs_1314</th>\n",
              "      <th>tma_facs_2172</th>\n",
              "      <th>tma_facs_1625</th>\n",
              "      <th>tma_facs_12420</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>Cts7</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Sdr42e1</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>46</td>\n",
              "      <td>52</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Commd5</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>89</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>50</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Fam170a</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Olfr748</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>...</th>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Bcar3</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>124</td>\n",
              "      <td>34</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>8</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Hk1</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>5</td>\n",
              "      <td>0</td>\n",
              "      <td>2</td>\n",
              "      <td>0</td>\n",
              "      <td>71</td>\n",
              "      <td>149</td>\n",
              "      <td>37</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>10</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>149</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>3</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Sh3glb1</th>\n",
              "      <td>18</td>\n",
              "      <td>0</td>\n",
              "      <td>10</td>\n",
              "      <td>0</td>\n",
              "      <td>6</td>\n",
              "      <td>678</td>\n",
              "      <td>205</td>\n",
              "      <td>79</td>\n",
              "      <td>86</td>\n",
              "      <td>4</td>\n",
              "      <td>...</td>\n",
              "      <td>187</td>\n",
              "      <td>0</td>\n",
              "      <td>408</td>\n",
              "      <td>69</td>\n",
              "      <td>12</td>\n",
              "      <td>0</td>\n",
              "      <td>11</td>\n",
              "      <td>0</td>\n",
              "      <td>13</td>\n",
              "      <td>36</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Calr4</th>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2410091C18Rik</th>\n",
              "      <td>10</td>\n",
              "      <td>0</td>\n",
              "      <td>118</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>130</td>\n",
              "      <td>36</td>\n",
              "      <td>1</td>\n",
              "      <td>0</td>\n",
              "      <td>...</td>\n",
              "      <td>182</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>24</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>2000 rows × 4000 columns</p>\n",
              "</div>"
            ],
            "text/plain": [
              "               tma_facs_14301  tma_facs_4756  tma_facs_14382  tma_facs_423  \\\n",
              "Cts7                        0              0               0             0   \n",
              "Sdr42e1                     0              0               0             0   \n",
              "Commd5                      0              0               0             0   \n",
              "Fam170a                     0              0               0             0   \n",
              "Olfr748                     0              0               0             0   \n",
              "...                       ...            ...             ...           ...   \n",
              "Bcar3                       0              0               0             0   \n",
              "Hk1                         0              0               5             0   \n",
              "Sh3glb1                    18              0              10             0   \n",
              "Calr4                       0              0               0             0   \n",
              "2410091C18Rik              10              0             118             0   \n",
              "\n",
              "               tma_facs_15612  tma_facs_4729  tma_facs_8174  tma_facs_9855  \\\n",
              "Cts7                        0              0              0              0   \n",
              "Sdr42e1                     0              0              0             46   \n",
              "Commd5                      0              0              0             89   \n",
              "Fam170a                     0              0              0              0   \n",
              "Olfr748                     0              0              0              0   \n",
              "...                       ...            ...            ...            ...   \n",
              "Bcar3                       0              0            124             34   \n",
              "Hk1                         2              0             71            149   \n",
              "Sh3glb1                     6            678            205             79   \n",
              "Calr4                       0              0              0              0   \n",
              "2410091C18Rik               0              0            130             36   \n",
              "\n",
              "               tma_facs_9055  tma_facs_8668  ...  tma_facs_5182  \\\n",
              "Cts7                       0              0  ...              0   \n",
              "Sdr42e1                   52              0  ...              0   \n",
              "Commd5                     0              0  ...              0   \n",
              "Fam170a                    0              0  ...              0   \n",
              "Olfr748                    0              0  ...              0   \n",
              "...                      ...            ...  ...            ...   \n",
              "Bcar3                      0              0  ...              0   \n",
              "Hk1                       37              0  ...              0   \n",
              "Sh3glb1                   86              4  ...            187   \n",
              "Calr4                      0              0  ...              0   \n",
              "2410091C18Rik              1              0  ...            182   \n",
              "\n",
              "               tma_facs_19106  tma_facs_11386  tma_facs_14006  tma_facs_7477  \\\n",
              "Cts7                        0               0               0              0   \n",
              "Sdr42e1                     0               0               0              0   \n",
              "Commd5                      0               0               0              0   \n",
              "Fam170a                     0               0               0              0   \n",
              "Olfr748                     0               0               0              0   \n",
              "...                       ...             ...             ...            ...   \n",
              "Bcar3                       0               0               0              0   \n",
              "Hk1                         0              10               0              0   \n",
              "Sh3glb1                     0             408              69             12   \n",
              "Calr4                       0               0               0              0   \n",
              "2410091C18Rik               0               0               0             24   \n",
              "\n",
              "               tma_facs_15352  tma_facs_1314  tma_facs_2172  tma_facs_1625  \\\n",
              "Cts7                        0              0              0              0   \n",
              "Sdr42e1                     0              0              0              0   \n",
              "Commd5                      0              0              0             50   \n",
              "Fam170a                     0              0              0              0   \n",
              "Olfr748                     0              0              0              0   \n",
              "...                       ...            ...            ...            ...   \n",
              "Bcar3                       0              8              0              0   \n",
              "Hk1                       149              0              0              3   \n",
              "Sh3glb1                     0             11              0             13   \n",
              "Calr4                       0              0              0              0   \n",
              "2410091C18Rik               0              0              0              0   \n",
              "\n",
              "               tma_facs_12420  \n",
              "Cts7                        0  \n",
              "Sdr42e1                     0  \n",
              "Commd5                      0  \n",
              "Fam170a                     0  \n",
              "Olfr748                     0  \n",
              "...                       ...  \n",
              "Bcar3                       0  \n",
              "Hk1                         0  \n",
              "Sh3glb1                    36  \n",
              "Calr4                       0  \n",
              "2410091C18Rik               0  \n",
              "\n",
              "[2000 rows x 4000 columns]"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "dataset"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "ea915110",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>0</th>\n",
              "      <th>1</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>tma_facs_14301</td>\n",
              "      <td>Monocyte</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>tma_facs_4756</td>\n",
              "      <td>Endothelial cell</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>tma_facs_14382</td>\n",
              "      <td>Monocyte</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>tma_facs_423</td>\n",
              "      <td>B cell</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>tma_facs_15612</td>\n",
              "      <td>Monocyte</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>...</th>\n",
              "      <td>...</td>\n",
              "      <td>...</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3995</th>\n",
              "      <td>tma_facs_15352</td>\n",
              "      <td>Monocyte</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3996</th>\n",
              "      <td>tma_facs_1314</td>\n",
              "      <td>B cell</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3997</th>\n",
              "      <td>tma_facs_2172</td>\n",
              "      <td>B cell</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3998</th>\n",
              "      <td>tma_facs_1625</td>\n",
              "      <td>B cell</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3999</th>\n",
              "      <td>tma_facs_12420</td>\n",
              "      <td>Epidermis</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "<p>4000 rows × 2 columns</p>\n",
              "</div>"
            ],
            "text/plain": [
              "                   0                 1\n",
              "0     tma_facs_14301          Monocyte\n",
              "1      tma_facs_4756  Endothelial cell\n",
              "2     tma_facs_14382          Monocyte\n",
              "3       tma_facs_423            B cell\n",
              "4     tma_facs_15612          Monocyte\n",
              "...              ...               ...\n",
              "3995  tma_facs_15352          Monocyte\n",
              "3996   tma_facs_1314            B cell\n",
              "3997   tma_facs_2172            B cell\n",
              "3998   tma_facs_1625            B cell\n",
              "3999  tma_facs_12420         Epidermis\n",
              "\n",
              "[4000 rows x 2 columns]"
            ]
          },
          "execution_count": 5,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "labels # Columns: cell ID and cell type"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "896fItX4vctP",
      "metadata": {
        "id": "896fItX4vctP"
      },
      "outputs": [],
      "source": [
        "n_types = 0 # To store the number of cell types present in the dataset, which is need to define the output layer of the model later on.\n",
        "\n",
        "# Convertes cell type labels from string ('B Cell', 'T Cell' etc) to integers (1, 2 etc)\n",
        "def convert_type2label(types):\n",
        "    global n_types\n",
        "    types = list(types)\n",
        "    unique_types = sorted(set(types)) # Sorting ensures consistent label order\n",
        "    n_types = len(unique_types)\n",
        "\n",
        "    type_to_label_dict = {t: i for i, t in enumerate(unique_types)}\n",
        "    labels = np.array([type_to_label_dict[t] for t in types])\n",
        "\n",
        "    return labels\n",
        "\n",
        "labels_int = convert_type2label(labels[1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "lld5Gr1ivgGg",
      "metadata": {
        "id": "lld5Gr1ivgGg"
      },
      "outputs": [],
      "source": [
        "# Reshaping labels from (n,) to (n, 1) to be compatible with DeepChem's Dataset class\n",
        "labels_int = labels_int.reshape(labels_int.shape[0],1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "RQCvRaSNQ7gC",
      "metadata": {
        "id": "RQCvRaSNQ7gC"
      },
      "outputs": [],
      "source": [
        "# gene IDs are stored\n",
        "train_genes = dataset.index.to_numpy()\n",
        "\n",
        "# genes x cells (features x samples) -> cells x genes (samples x features)\n",
        "dataset = np.transpose(dataset)\n",
        "\n",
        "# Creating Deepchem 'Dataset' object\n",
        "dataset = dc.data.NumpyDataset(X=np.array(dataset,dtype=np.float32),y = labels_int, ids= labels[0])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "B0TK7Uz-vxmp",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 304
        },
        "id": "B0TK7Uz-vxmp",
        "outputId": "cfdcdf5f-e0ad-46c0-d2b7-441c496045c5"
      },
      "outputs": [],
      "source": [
        "splitter = dc.splits.SingletaskStratifiedSplitter()\n",
        "train_dataset, test_dataset = splitter.train_test_split(dataset)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2a5a810e",
      "metadata": {
        "id": "2a5a810e"
      },
      "source": [
        "### 2. Different sources"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b2509e46",
      "metadata": {
        "id": "b2509e46"
      },
      "source": [
        "When using datasets from different sources, the gene sets may not completely overlap. Since the ACTINN model defines its layers based on the gene set, it's essential that both the training and testing sets contain the same genes. To ensure this, we first identified the genes common to both datasets and filtered out those that were not shared.\n",
        "\n",
        "We will now demonstrate how to prepare data from different sources, which is the method we'll use for the rest of this tutorial"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "cd1aa909",
      "metadata": {
        "id": "cd1aa909"
      },
      "outputs": [],
      "source": [
        "train_set = pd.read_hdf(os.path.join(current_dir,'train_sample.h5'))\n",
        "train_labels = pd.read_csv(os.path.join(current_dir,'train_label_sample.txt'), header=None, sep='\\t')\n",
        "test_set = pd.read_hdf(os.path.join(current_dir,'test_sample.h5'))\n",
        "test_labels = pd.read_csv(os.path.join(current_dir,'test_label_sample.txt'), header=None, sep='\\t')\n",
        "\n",
        "common_genes = train_set.index.intersection(test_set.index)\n",
        "common_genes = sorted(common_genes)\n",
        "\n",
        "train_set = train_set.loc[common_genes]\n",
        "test_set = test_set.loc[common_genes]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5559ab37",
      "metadata": {
        "id": "5559ab37"
      },
      "source": [
        "Convertes cell type labels from string ('B Cell', 'T Cell' etc) to integers (1, 2 etc)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "4cbe0639",
      "metadata": {
        "id": "4cbe0639"
      },
      "outputs": [],
      "source": [
        "train_labels_int = convert_type2label(train_labels[1])\n",
        "test_labels_int = convert_type2label(test_labels[1])\n",
        "\n",
        "# Reshaping labels from (n,) to (n, 1) to be compatible with DeepChem's Dataset class\n",
        "train_labels_int = train_labels_int.reshape(train_labels_int.shape[0],1)\n",
        "test_labels_int = test_labels_int.reshape(test_labels_int.shape[0],1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "cZjtIemOzZAR",
      "metadata": {
        "id": "cZjtIemOzZAR"
      },
      "outputs": [],
      "source": [
        "# gene IDs are stored\n",
        "train_genes = train_set.index.to_numpy()\n",
        "test_genes = test_set.index.to_numpy()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "28VzJlRT5JhT",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "28VzJlRT5JhT",
        "outputId": "0a7e7402-fdb9-46c1-ab82-f17323fac9ad"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "2000 2000\n"
          ]
        }
      ],
      "source": [
        "print(len(train_genes),len(test_genes))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "J39tJUqYzFyq",
      "metadata": {
        "id": "J39tJUqYzFyq"
      },
      "outputs": [],
      "source": [
        "# genes x cells (features x samples) -> cells x genes (samples x features)\n",
        "train_set = np.transpose(train_set)\n",
        "test_set = np.transpose(test_set)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "X3AFcLECylTJ",
      "metadata": {
        "id": "X3AFcLECylTJ"
      },
      "outputs": [],
      "source": [
        "# Creating Deepchem 'Dataset' object\n",
        "train_dataset = dc.data.NumpyDataset(X=np.array(train_set,dtype=np.float32),y = train_labels_int, ids= train_labels[0])\n",
        "test_dataset = dc.data.NumpyDataset(X=np.array(test_set,dtype=np.float32), y = test_labels_int, ids= test_labels[0])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "W9qpsgCn5n8I",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "W9qpsgCn5n8I",
        "outputId": "184a79ca-af50-4a54-9731-728bc92244c4"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<NumpyDataset X.shape: (4000, 2000), y.shape: (4000, 1), w.shape: (4000, 1), task_names: [0]>"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "train_dataset"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "rHDCAgrx5plj",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "rHDCAgrx5plj",
        "outputId": "f4daf5ce-cf64-469c-909f-7be84f9d767d"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<NumpyDataset X.shape: (1000, 2000), y.shape: (1000, 1), w.shape: (1000, 1), task_names: [0]>"
            ]
          },
          "execution_count": 17,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "test_dataset"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "272a3e17",
      "metadata": {
        "id": "272a3e17"
      },
      "source": [
        "## Data Preprocessing\n",
        "\n",
        "Starting from raw counts, a scRNA-Seq data analysis typically includes normalization, feature selection, and dimension reduction steps.\n",
        "\n",
        "In a standard RNA-seq workflow, RNA molecules are first captured from cells, reverse transcribed(i.e RNA is converted to DNA) into complementary DNA and sequenced. The resulting short reads are then computationally aligned to reference genes to obtain count data. Each stage of this process introduces some degree of technical variability, even among cells that are biologically identical.\n",
        "\n",
        "The preprocessing pipeline includes the following steps:\n",
        "1) Library-size normalize to 10,000 counts per cell and Log2(x+1) transform\n",
        "2) Filter genes by total expression (1st–99th percentile)\n",
        "3) Filter genes by coefficient of variation (1st–99th percentile)\n",
        "\n",
        "**Normalization** seeks to adjust for differences in experimental conditions between samples (individual cells), so that these do not confound true biological differences.[6]\n",
        "This step aims to adjust the raw counts in the dataset for variable sampling effects by scaling the observable variance to a specified range. Several normalization techniques are used in practice varying in complexity. For our analysis, we'll be using the shifted normalisation method as described in [7]. Log transformation reduces data skewness, improving suitability for machine learning applications.\n",
        "\n",
        "The image below illustrates how the distribution of counts changes after applying the shifted logarithm, compared to the total counts in the raw dataset\n",
        "\n",
        "<img src=\"https://www.sc-best-practices.org/_images/e7db84d20620d812e8d3b77a196d247d0ac339a306ac38eacb69c4a0a69f0321.png\">\n",
        "\n",
        "**Feature selection**, or identification of informative genes, is accomplished by ranking genes using total expression and coefficient of variation(std/mean), followed by the removal of the top and bottom 1%. Highly expressed genes might be housekeeping genes(genes that are essential for basic cellular functions and expressed in all cells). Genes that have low counts are likely to be technical noise. Genes with a low CV have very stable expression levels across all cells relative to their mean, like housekeeping genes, their lack of variation makes them uninformative. Extremely high CV values are often artifacts associated with genes that have a very low average expression. For these genes, the detection of just a few transcripts in only a handful of cells can create a large standard deviation relative to a near-zero mean, resulting in an artificially inflated CV\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c8e5c802",
      "metadata": {},
      "source": [
        "## Training Data Preprocessing"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "df28cbc0",
      "metadata": {},
      "source": [
        "library-size normalize to 10,000 and log2(x + 1) transform"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "48862c97",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(4000, 2000)"
            ]
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# extract gene names & data array\n",
        "gene_names = train_genes\n",
        "X = train_dataset.X\n",
        "\n",
        "row_sums = X.sum(axis=1, keepdims=True)  # shape (n_cells, 1)\n",
        "X /= row_sums # broadcast divide\n",
        "X *= 10000\n",
        "\n",
        "X.shape # (number of cells, number of genes)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "662c879c",
      "metadata": {},
      "source": [
        "The counts are increased by 1, and the log2 value is calculated"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "f1eeb4c2",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "f1eeb4c2",
        "outputId": "78159223-b246-4721-e0f5-3bba97c5c8ee"
      },
      "outputs": [],
      "source": [
        "# log2(x + 1) transform\n",
        "X = np.log2(X + 1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5ff8d469",
      "metadata": {},
      "source": [
        "To filter out outlier genes, the genes with the highest 1% and lowest 1% expression were removed."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "9b16db0b",
      "metadata": {
        "id": "9b16db0b"
      },
      "outputs": [],
      "source": [
        "expr = X.sum(axis=0)\n",
        "low, high = np.percentile(expr, [1, 99])\n",
        "mask_expr = (expr >= low) & (expr <= high)\n",
        "\n",
        "X = X[:, mask_expr]\n",
        "gene_names = gene_names[mask_expr]\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "973536f9",
      "metadata": {},
      "source": [
        "To remove genes with zero mean (genes that aren't expressed in any of the cells)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "c958abbb",
      "metadata": {},
      "outputs": [],
      "source": [
        "mean_expr = X.mean(axis=0)\n",
        "mask_mean = mean_expr > 0\n",
        "X = X[:, mask_mean]\n",
        "gene_names = gene_names[mask_mean]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "87378bf8",
      "metadata": {
        "id": "87378bf8"
      },
      "source": [
        "The genes with the highest 1% and the lowest 1% standard deviation are removed."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "cf0e30bd",
      "metadata": {
        "id": "cf0e30bd"
      },
      "outputs": [],
      "source": [
        "mean_expr = X.mean(axis=0)\n",
        "cv = X.std(axis=0) / mean_expr\n",
        "low_cv, high_cv = np.percentile(cv, [1, 99])\n",
        "mask_cv = (cv >= low_cv) & (cv <= high_cv)\n",
        "\n",
        "X = X[:, mask_cv]\n",
        "gene_names = gene_names[mask_cv]\n",
        "train_genes = gene_names"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "id": "sevCNmAo9TLh",
      "metadata": {
        "id": "sevCNmAo9TLh"
      },
      "outputs": [],
      "source": [
        "# Creating Deepchem 'Dataset' object\n",
        "train_dataset = dc.data.NumpyDataset(X=X,y = train_labels_int, ids= train_labels[0])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "7b480c75",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<NumpyDataset X.shape: (4000, 1865), y.shape: (4000, 1), w.shape: (4000, 1), task_names: [0]>"
            ]
          },
          "execution_count": 24,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "train_dataset"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "471975c9",
      "metadata": {
        "id": "471975c9"
      },
      "source": [
        "## Test Data Preprocessing\n",
        "\n",
        "The exact same normalisation aproach used for train set is used for test set.\n",
        "\n",
        "First, each cell’s expression value was normalized to its total expression value and multiplied by a scale factor\n",
        "of 10 000. The counts were increased by 1, and the log2 value is calculated. \n",
        "\n",
        "The gene list obtained from the train set filtering steps is used to mask the test set. In the original ACTINN implementation, genes were filtered using both the training and test sets. To avoid potential information leakage, we instead apply filtering based solely on the training set."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "id": "ac7f8796",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ac7f8796",
        "outputId": "fe51a78d-2aea-457a-b6fa-2e074ab466e1"
      },
      "outputs": [],
      "source": [
        "# extract data array\n",
        "X = test_dataset.X\n",
        "\n",
        "# library-size normalize to 10,000 (in-place)\n",
        "row_sums = X.sum(axis=1, keepdims=True)  # shape (1, n_cells)\n",
        "X /= row_sums  # broadcast divide\n",
        "X *= 10000\n",
        "\n",
        "# log2(x + 1) transform\n",
        "X = np.log2(X + 1, out=X)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "id": "prjpO0pG35wU",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "prjpO0pG35wU",
        "outputId": "a1bab2ea-4b4e-4d1e-ef07-0254b71e0a2a"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1865 2000\n"
          ]
        }
      ],
      "source": [
        "print(len(train_genes), len(test_genes))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "07ccf58c",
      "metadata": {},
      "source": [
        "The gene list obtained from the train set filtering steps is used to mask the test set."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "27ebb3ed",
      "metadata": {
        "id": "27ebb3ed"
      },
      "outputs": [],
      "source": [
        "test_gene_mask=[i in train_genes for i in test_genes]\n",
        "X = X[:,test_gene_mask]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "id": "8GfNQ-tP7xyu",
      "metadata": {
        "id": "8GfNQ-tP7xyu"
      },
      "outputs": [],
      "source": [
        "test_dataset = dc.data.NumpyDataset(X=np.array(X,dtype=np.float32), y=test_labels_int, ids=test_labels[0])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "id": "KFy_Eka0-fUh",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KFy_Eka0-fUh",
        "outputId": "88d3b2cf-61b9-4d6f-a1ac-896ce874b249"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<NumpyDataset X.shape: (1000, 1865), y.shape: (1000, 1), w.shape: (1000, 1), task_names: [0]>"
            ]
          },
          "execution_count": 29,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "test_dataset"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0b2bbf47",
      "metadata": {
        "id": "0b2bbf47"
      },
      "source": [
        "# Model definition\n",
        "\n",
        "ACTINN is implemented with deepchem Model class as wrapper. It uses a neural network that containing an input layer, three hidden layers and an output layer. The input layer has a number of nodes equal to the number of genes in the training set. The 3 hidden layers have 100, 50 and 25 nodes, respectively. The output layer has a number of nodes equal to the number of cell types in the training set. For the activation functions, we used the softmax function for the ouput layer and the rectified linear unit (ReLU) function for the other layers. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e5124361",
      "metadata": {
        "id": "e5124361"
      },
      "outputs": [],
      "source": [
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "from deepchem.models.torch_models import TorchModel\n",
        "from deepchem.models.losses import SparseSoftmaxCrossEntropy\n",
        "from deepchem.models.optimizers import Adam\n",
        "from deepchem.models.optimizers import ExponentialDecay\n",
        "from deepchem.metrics import from_one_hot\n",
        "\n",
        "\n",
        "class ActinnClassifier(nn.Module):\n",
        "\n",
        "    def __init__(self, output_dim=None, input_size=None):\n",
        "\n",
        "        if output_dim == None or input_size == None:\n",
        "            raise ValueError('Must explicitly declare input dim (num features) and output dim (number of classes)')\n",
        "\n",
        "        super(ActinnClassifier, self).__init__()\n",
        "        self.inp_dim = input_size\n",
        "        self.out_dim = output_dim\n",
        "\n",
        "        # feed forward layers\n",
        "        self.classifier = nn.Sequential(\n",
        "                                        nn.Linear(self.inp_dim, 100),\n",
        "                                        nn.ReLU(),\n",
        "\n",
        "                                        nn.Linear(100, 50),\n",
        "                                        nn.ReLU(),\n",
        "\n",
        "                                        nn.Linear(50, 25),\n",
        "                                        nn.ReLU(),\n",
        "\n",
        "                                        nn.Linear(25, output_dim)\n",
        "                                        )\n",
        "\n",
        "    def forward(self, x):\n",
        "        \"\"\"\n",
        "        Forward pass of the classifier\n",
        "        \"\"\"\n",
        "        logits = self.classifier(x)\n",
        "        probabilities = F.softmax(logits, dim=1)\n",
        "        predictions = from_one_hot(probabilities.cpu().detach(), axis=1)\n",
        "        return (predictions, logits)\n",
        "    \n",
        "\n",
        "class ACTINNModel(TorchModel):\n",
        "    def __init__(self, output_dim = None, input_size = None, **kwargs):\n",
        "\n",
        "        self.model = ActinnClassifier(output_dim, input_size)\n",
        "        self.criterion = SparseSoftmaxCrossEntropy()\n",
        "        cf_optimizer = Adam(learning_rate=0.0001,\n",
        "                            beta1=0.9,\n",
        "                            beta2=0.999,\n",
        "                            epsilon=1e-08,\n",
        "                            weight_decay=0.005,\n",
        "                            )\n",
        "\n",
        "        cf_decayRate = 0.95\n",
        "        cf_lr_scheduler = ExponentialDecay(initial_rate=0.0001, decay_rate=cf_decayRate, decay_steps=1000)\n",
        "        super(ACTINNModel,\n",
        "              self).__init__(self.model,\n",
        "                             loss=self.criterion,\n",
        "                             optimizer=cf_optimizer,\n",
        "                             learning_rate=cf_lr_scheduler,\n",
        "                             output_types=['prediction', 'loss'],\n",
        "                             **kwargs)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "88050735",
      "metadata": {
        "id": "88050735"
      },
      "source": [
        "## Model training and evaluation"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "id": "815b93bf",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "815b93bf",
        "outputId": "349ccde9-5306-4051-e4aa-eb095927dcc8"
      },
      "outputs": [],
      "source": [
        "model = ACTINNModel(output_dim= n_types,input_size= len(train_genes))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 161,
      "id": "5ad7fa1b",
      "metadata": {
        "id": "5ad7fa1b"
      },
      "outputs": [],
      "source": [
        "model.fit(train_dataset, nb_epoch=7)\n",
        "predictions = model.predict(test_dataset)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 162,
      "id": "58020421",
      "metadata": {},
      "outputs": [],
      "source": [
        "classification_metric = dc.metrics.Metric(dc.metrics.accuracy_score)\n",
        "scores = model.evaluate(test_dataset, [classification_metric], n_classes=n_types)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 163,
      "id": "595c4d0a",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "0.971\n"
          ]
        }
      ],
      "source": [
        "print(scores['accuracy_score'])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0990cfb6",
      "metadata": {
        "id": "0990cfb6"
      },
      "source": [
        "### References\n",
        "[1] Tang, F., Barbacioru, C., Wang, Y., Nordman, E., Lee, C., Xu, N., ... & Surani, M. A. (2009). mRNA-Seq whole-transcriptome analysis of a single cell. Nature methods, 6(5), 377-382.\n",
        "\n",
        "[2] Alberts, B., Johnson, A., Lewis, J., Morgan, D., Raff, M., Roberts, K., & Walter, P. (2014). Molecular Biology of the Cell (6th ed.). Garland Science.\n",
        "\n",
        "[3] Jorge A. Tzec-Interián, Daianna González-Padilla, Elsa B. Góngora-Castillo,.(2025). Bioinformatics perspectives on transcriptomics: A comprehensive review of bulk and single-cell RNA sequencing analyses.\n",
        "\n",
        "[4] Feiyang Ma , Matteo Pellegrini. (2019), ACTINN: automated identification of cell types in single cell RNA sequencing. Bioinformatics Volume 36, Issue 2, January 2020, Pages 533–538.\n",
        "\n",
        "[5] Hwang, B., Lee, J. H., & Bang, D. (2018). Single-cell RNA sequencing technologies and bioinformatics pipelines. Experimental & Molecular Medicine, 50(8), 96. https://doi.org/10.1038/s12276-018-0071-8\n",
        "\n",
        "[6] Townes, F. W., Hicks, S. C., Aryee, M. J., & Irizarry, R. A. (2019). Feature selection and dimension reduction for single‑cell RNA‑Seq based on a multinomial model. Genome Biology, 20(1), Article 295. https://doi.org/10.1186/s13059-019-1861-6 \n",
        "\n",
        "[7] Ahlmann‑Eltze, Constantin, and Wolfgang Huber. 2023. “Comparison of Transformations for Single‑Cell RNA‑Seq Data.” Nature Methods 20 (5): 665–672. https://doi.org/10.1038/s41592‑023‑01814‑1. \n",
        "\n",
        "[8] https://figshare.com/articles/dataset/ACTINN/8967116\n",
        "\n",
        "[9] Lambrechts, Diether, et al. Single-cell transcriptomics of human and mouse lung cancers reveals conserved myeloid populations across individuals and species. Nature Biotechnology 36.10 (2018): 1–13. https://doi.org/10.1038/nbt.4220\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2b0df63d",
      "metadata": {},
      "source": [
        "# Congratulations! Time to join the Community!\n",
        "\n",
        "Congratulations on completing this tutorial notebook! If you enjoyed working through the tutorial, and want to continue working with DeepChem, we encourage you to finish the rest of the tutorials in this series. You can also help the DeepChem community in the following ways:\n",
        "\n",
        "## Star DeepChem on [GitHub](https://github.com/deepchem/deepchem)\n",
        "This helps build awareness of the DeepChem project and the tools for open source drug discovery that we're trying to build.\n",
        "\n",
        "## Join the DeepChem Discord\n",
        "The DeepChem [Discord](https://discord.gg/cGzwCdrUqS) hosts a number of scientists, developers, and enthusiasts interested in deep learning for the life sciences. Join the conversation!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b121bf89",
      "metadata": {},
      "outputs": [],
      "source": [
        "@manual {Bioinformatics,\n",
        " title={Cell type identification in scRNA-seq data using ACTINN and Deepchem},\n",
        " organization={DeepChem},\n",
        " author={Adithya, Harindhar and Singh, Rakshit kr. and  Ramsundar, Bharath},\n",
        " howpublished = {\\url{https://github.com/deepchem/deepchem/blob/master/examples/tutorials/Cell_type_identification_using_scRNAseq_data.ipynb}},\n",
        " year={2025},\n",
        "}"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "dc_new",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.8"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}
